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Abstract 

This paper presents a numerical study of three-dimensional laminar mixed convec- 
tion within a liquid flowing on a horizontal channel heated uniformly from below. 
The upper surface is free and assumed to be flat. The coupled Navier-Stokes and 
energy equations are solved numerically by the finite volume method taking into 
account the thermocapillary effects (Marangoni effect). When the strength of the 
buoyancy, thermocapillary effects and forced convective currents are comparable 



1 Corresponding author: kamal.elomari@univ-pau.fr 
Preprint submited to International Journal of Heat and Fluid Flow 



1 



1 Introduction 



2 



(Ri ^2 0(1) and Bd = Ra/Ma ±2 0(1)), the results show that the development of 
instabilities in the form of steady longitudinal convective rolls is similar to those 
encountered in the Poiseuille-Rayleigh-Benard flow. The number and spatial dis- 
tribution of these rolls along the channel depend on the flow conditions. The ob- 
jective of this work is to study the influence of parameters, such as the Reynolds, 
Rayleigh and Biot numbers, on the flow patterns and heat transfer characteris- 
tics. The effects of variations in the surface tension with temperature gradients 
(Marangoni effect) are also considered. 

Keywords: mixed convection ; longitudinal rolls ; Marangoni effect ; 
horizontal channel ; free surface. 



1 Introduction 

The study of mixed convection in a horizontal channel with a free upper 
surface is of considerable interest for many industrial applications, such as 
heat exchangers, evaporative cooling devices, and chemical vapor process. In 
particular conditions, the resulting buoyancy forces induce the development 
of longitudinal vortex rolls, which constitute the first stage of transition to 



turbulence and contribute to an enhancement of heat transfer f lmura et al. 



19781 ) 



Several analytical, numerical and experimental analyses involving mixed con- 
vection heat transfer in different channel geometries have been reported in 
the literature. However, it is common practice to model a thin liquid film 
on a wetted wall as a boundary condition for the heat and mass transfer 
in the adjacent gas flow. Thus, laminar mixed conv ection hea t and mass 



transfer in wetted ve rtical ducts has been analyzed by iLin et al.l ( 119881 ) and 



Tsay and Yanl (1199 ll ) . They showed that the combined thermal and mass dif- 



fusion buoyancy forces have a great effect for sm all Reynolds number flows 



with high inlet temperatures. lOulaid et al.l (120 lOh conducted a similar study 



that analyzed the effect of flow reversal on inducing flow instability. When 
the thickness of the liquid film is not negligible, it became interesting to 
study the flow mechanisms inside the film itself. 

In this study, we consider the case of mixed convection inside a liquid 
film with a free surface where thermocapillary effects are present. When a 
free liquid surface is present, the surface tension variations resulting from the 
temperature gradient along the surface may induce motion within the fluid, 
called thermocapillary flow (i.e., thermal Marangoni convection). To the best 
of our knowledge, the combination of buoyancy, thermocapillary forces and 
forced axial flow has not been addressed in the literature. However, this 
original flow situation is the combination of other flow situations (of physical 
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effects) that has been extensively studied. One of these situations is the flow 
inside a horizontal channel with the superposition of a forced axial flow and 
vertical natural convection currents. This unstable flow situation is known as 
the Poiseuille-Rayleigh-Benard (PRB) flow and has flow structures that are 
similar to the flow along an open channel that we are interested in studying. 

Many studies, both theoretical and experimental, have been conducted on 
mixed convection heat transfer in horizontal channels with rigid upper sur- 
face. The onset of a secondary flow, in the form of lo ngitudinal rolls, takes 
place around a critical Ra yleigh number of Ra c = 1708 (lAkivama et all 1 1971 : 



Yasuo and Yutakal. 119661) for which the r oll dia meter is nearly equal to the 



channel height. lOstrach and Kamotanil (119751 ) performed experiments for 



fully developed air flow between isothermal plates and noted that the vortex 

" feoooh 



Nicolas et al. 



rolls become irregular for Re = 38 and Ra > 8000. 
presented the linear stability analysis of a ful ly developed Poiseuill e flow in 
a duct and established a stability diagram. iBonnefoi et al.l ( 120041 1 investi- 
gated the thermoconvective instabilities for a mixed convection water flow 
in a horizontal rectangular duct uniformly heated from below. The authors 
distinguish several regim es of flow a n d, sp ecifically, the domains of ther- 
moconvective instability. IWang et al.l ( 120031 ) investigated three-dimensional 
mixed convection flow and heat transfer in a horizontal duct using constant 
or variable thermophysical property models. Th e numerical results from the 
two different approaches have been compared. iLuijkx et al.l ( I198ll ) demon- 



strated the existence of unsteady transverse thermoconvective rolls with the 
axis perpendicular to the main flow for the flow of silicone oil (Pr = 450). 
For low Reynolds numbers, the critical Rayleigh number corresponding to 
the onset of the transverse rolls was found to be a function of the aspect 
ratio and the Prandtl number. A transition from transverse to longitudinal 
roll flow structure was found, and the corresp onding boundar y in the (Re - 



Ra) plane was exp erimentally es t ablish ed by lOuazzani et al.l (Il990l Il993l ) 



and numerically by iNicolas et al.l ( 119971 ) . A comprehensive literature review 
dedicated to t he Poiseuille-Rayleigh-Benard flow structure was presented by 



Nicolasl (120021 1 . 



A second flow situation that is similar to the one studied in this paper 
is the case of a channel whose upper surface is free but without thermocap- 
ill ary forces. This cas e has been less studied. Experiments were performed 
by iGilpin et al.l (119781 ) to confirm the occurrence and growth of longitudinal 
vortices in the laminar boundary layer developing in the water over a heated 
horizontal flat plate with a uniform surface temperature. Temperature pro- 
files across the boundary layer were measured for flows with and without 
vortices to show qualitatively the effect that longitudinal vortices have on 
the heat transfer rate on the plate. An experimental study was also per- 
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formed by IWang et al.l (119831 ) to determine the hydrodynamic and thermal 
conditions of water flow in an open channel that is uniformly heated from 
below. The authors also used a two-dimensional boundary layer model to 
predict the flow behavior. 

Flow situation studies, where thermocapillary forces (Marangoni convec- 
tion) have been taken into account, are essentially related to flows without 
forced convection. These flows have been extensive l y stud i ed in the literature 



f|Goncharova an d Kabov, 201 01; ISimanovskii et all 120081 : ISadiq et all 12010 



Mao et al.l . 120081 : iMedale and Cochelinl . 120091). A pionee r study of thermo- 

Pearsonl (119581 ) who analyzed the 



capillary instabilities was conducted by 
stability of a horizontal layer with a non-deformable free surface. The liq- 
uid phase is modeled using the classical fluid mechanics equations, while the 
gas phase is included in the model through a boundary condition at the 
liquid-gas interface, which is characterized by a Biot number. Through a lin- 
ear stability analysis, Pearson identified the existence of a critical threshold 
and determined a critical Marangoni number from which the liquid layer be- 
comes unstable due to the thermocap il lary c onvective motion that appears 



spontaneously. IScriven and Sternlingi (119641 ) extended Pearson's study by 



considering a deformable liquid-gas interface. Unlike the situation addressed 
by Pearson, they concluded that the liquid layer is still unstable and there 
is no critical Marangoni number. More recently, the coupling between ther- 
mocapillary and buoyancy has been consi dered through different a pproaches. 
To verify the influence of buoyancy forces, lLaure and Rou x ( 1989 ) conducted 



severa l studies for low Prandtl number ( Pr) fluids, while Parmentier et al. 
(Il993[ ) studied liquids with Pr up to 10. IBurguete et al.l ( 120011 ) found that 
thermocapillary-buoyancy convection could destabilize the flow into differ- 
ent patterns, depending on the temperature difference AT and on the liq- 
uid pool depth d. For small d values, the system exhibits hydrothermal 
waves, while for larger d values, stationary longitudinal rolls are observed. 
The stability of buoyant-the rmocapillary instabilities has been analyzed by 
Mercier and Normandl (Il996| ). They introduce a Biot number to describe the 



heat transfer at the top free surface and showed the existence of a transition 
from oscillatory to stationary modes when the Bond number (the ratio be- 
tween the Rayleigh and the Marangoni numbers) is increased. This transition 
depends on the Biot number. 

Other studies have examined a thin liquid film flowing ove r an inclined 
heated plate, taking into account o nly the thermocapillary effect ( ISadiq and Ushal . 



2005t iThiele and Knoblochl . 120041 ) . In these studies, most authors did not 



consider the thermoconvective instabilities that occur within a liquid film 
uniformly heated from below. Indeed, when the strength of buoyancy, ther- 
mocapillary and forced convective currents are comparable, (Ri ^ 0(1) and 
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Bd = Raj Ma ±2 0(1)), and the instabilities develop in the form of longitu- 
dinal convective rolls similar to those found in Poiseuille-Rayleigh-Benard. 

This study is dedicated to the numerical simulation of three-dimensional 
laminar mixed convection of a liquid in a horizontal channel uniformly heated 
from below with a free upper surface. The objective is to study the influence 
of the control parameters (i.e., Reynolds, Rayleigh and Biot numbers) on the 
heat transfer in the liquid. The effects of the surface tension and their varia- 
tions with temperature gradients (Marangoni effect) are taken into account. 
The paper is organized as follows. The physical model and governing equa- 
tions are first described, followed by a description of the numerical resolution 
method and code validation. The effect of the studied parameters on the lon- 
gitudinal rolls and heat transfer are then assessed. Finally, the impact of the 
coupled effects of natural convection and Marangoni convection is analyzed. 
Conclusions are then drawn based on these results. 

2 Physical model 

We consider a liquid flow (Pr = 7) in a horizontal rectangular channel of 
height H, width 10H and length 50H with a free upper surface (Figured]). 
We define the transverse aspect ratio as T = width/height = 10. The side 
walls of the channel are adiabatic. To avoid the formation of the convective 
rolls immediately at the inlet section, we consider an adiabatic entrance over 
a length of 2H. Beyond this adiabatic entrance, the bottom-wall is kept at a 
constant and uniform temperature T c , which is higher than that of the inlet 
liquid T (T c > T ). The temperature of the ambient air over the liquid flow 
is fixed at T . 



Fig. 1: The studied configuration and corresponding boundary conditions. 

The fluid is assumed to be Newtonian and incompressible and under- 
goes steady laminar flow. The thermophysical properties are assumed to be 
constant except for the density in the buoyancy term and the surface ten- 
sion. The Boussinesq assumption is assumed, with p taken to be constant 



adiabatic side walls 




adiab. 
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everywhere, except in the buoyancy term, where a linear law is considered: 
p = p (l — f3(T — T )), where (3 is the thermal expansion coefficient and p 
and T are the reference values for the density and temperature, respectively. 
This assumptio n is justified for a suffici ently weak temperature difference 
across the fluid (IGrav and Giorginil . Il976l ) . The variation of the surface ten- 
sion with temperature is modeled by the following linear approximation: 
a = (j — j(T — T ), where a is the surface tension at temperature T and 
7 = —da/dT. The free surface is assumed to be flat and is subject to a 
convective heat exchange with ambient air characterized by a heat transfer 
coefficient h. The surface deformation can be neglected if the Crispation 
number Cr = pvajoR <C 1 (jDavisl . 119871 ). where p, v and a denote the fluid 
density, kinematic viscosity and thermal diffusivity, respectively. Moreover, 
the Galileo number Ga = gH 3 /ua 3> 1 characterizes the relative impor- 
tance of gravity and diffusion (g is the gravity constant). A large value of 
the Galileo number indicates that gravity stabilizes the long-wave mode of 
the free surface deformation. Such conditions are shown be satisfied for the 
problem considered here, where Cr = 1.42 x 10 -13 and Ga = 7 x 10 13 for 
Ra = 5000, Re = 15, Bi = 15 and Ma = 500. Using these assumptions, 
the equations governing the conservation of mass, momentum and energy are 
written as follows: 



V.v 



v.VT 







(v.V)v = Vp + uV 2 v + Pg(T -T ) 

Po 



aV 2 T 



(1) 
(2) 

(3) 



Where v denote the three-dimensional velocity field and T is the temperature 
field. 

At the bottom of the channel, the velocity satisfies the no-slip condition, 
and the wall is assumed to be adiabatic for x G [— 2H, 0[ and isothermal 
(T = T c ) for x G [0,48-ffj: 

dT 

v = 0, — = for x e [-2H, 0[; T = T c for x G [0, A8H] at z = 0. (4) 

The velocity satisfies also the no-slip condition for the side walls, which 
are assumed to be adiabatic: 

dT 

v = 0, — = at y = and y = IQH. (5) 
dy 
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By assuming a planar surface, a shear stress boundary condition is im- 
posed at the free surface, derived from the balance between the surface ten- 
sion forces and the viscous stresses in the fluid: 

dv x dadT dv v da dT 

The other boundary condition at the interface involves the vertical com- 
ponent of the velocity v z (non-permeability condition): 

Vz = o at z = H. (7) 
The thermal boundary condition at the free surface is: 
dT 

- X— = h(T - T ) at z = H. (8) 

oz 

The first term represents the heat conduction in the liquid, with A denot- 
ing the fluid thermal conductivity, while the second term expresses the heat 
flux density to the ambient air considered at To. 

Finally, at the inlet, we use a half Poiseuille profile for velocity. This 
choice can be justified by considering a slip condition at the upper surface 
(lUf = = 0). The three-dimensional Poiseuille profile in the case of 
a rectangular chann el with rigid walls has been computed analytically in 



Nicolas et al.l ( 120001 ). The inlet temperature is fixed at Tq: 



v x = vi PmsM ; Vy = v z = 0; T = T at x = -2H. (9) 

Using the channel height H, the mean flow velocity U m and pU^ as refer- 
ence quantities for the lengths, velocities and pressure, respectively, and using 
the reduced temperature 9 = (T — T )/(T C — T ), the governing equations 
take the following dimensionless form: 



v.y = o (io) 

1 Rn 

{ v.v)v = -vp + -Vv + — 2 e% (n) 

™" = pk f2e (12) 

Thus, the dimensionless parameters that appear in these equations are 
the Reynolds number Re = ^JL^ the Rayleigh number Ra = 9l3 ^ Tc ~ T ^ H an d 
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the Prandtl number Pr = -. 

a 

The boundary conditions (4) to (9) become: 
86 

y = 0, — = for X G [-2, 0[, = 1 for X G [0, 48] at Z = 0. (13) 

C/Z/ 

V = 0, qy =0 at F = and y = 10. (14) 

dV x Made dVy MadO 

-dZ = -T~edX' -dZ=-T~edY' Vz = ° Z = 1. (15) 

f)0 

— = -Bid at Z= 1. (16) 

= Vi PoU(YiZ) , V Y = V z = 0, 6 = at X = -2. (17) 

Where Ma = ^^ LK , Bi = and Pe = are the Marangoni, Biot and 
Peclet numbers, respectively. To quantify the heat transfer at the bottom 
wall, we define the local, the average and the transverse average Nusselt 
numbers by: 

N Ul {X,Y) = -(^\ , Nu M , = ~jNu l dS 

Z= ° s (18) 

10 



Nu tav {X) = ±- [ Nu t dY 

1U JO 



with S is the bottom heated surface of the channel. 



3 Numerical method 

The governing equations ([T], [2] and ED are solved using an in-house code 
called Tamaris. This code has a three-dimensional unstructured finite- volume 
framework that is applied to hybrid meshes. Variable values (v, p and T) are 
stored at cell centers in a collocated arrangement. Cell shapes can be of dif- 
ferent forms (tetrahedral, hexahedral, prismatic or pyramidal). To describe 
the discretization method used in the code, we can write Eqs.(J2]) and (J3J) in 
the generic convection-diffusion form with respect to a conserved variable 0: 

Jp<j)dV + J P (f)U.ndS = jAV(f).n,dS + Js^dV (19) 

V S S V 

Where A is a diffusion coefficient and is a source term. The spatial schemes 
for approximating the diffusive and convective fluxes are both second-order 
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accurate. The diffusion term is discretized by approximating the surface 
integrals with a sum over all cell faces /: 

[AV<p.ndS = J2 A f A f(V ( i ) )f ft f ( 20 ) 

s f 

where Af is the area of the face /. For unstructured meshes, orthogonality is 
an exception that needs to be handled correctly Thus, the normal gradient 
(y<p) f-fif is decomposed into an implicit contribution that uses the values 
of 4> at the centers of the two cells sharing face / (the first term on the RHS 
of Eq. (jZTJ)) and a non-orthogonality correction term treated explicitly by a 
deferred approach to preserve the second-order accuracy of the centered dif - 



(V0) r n/ = + V0. ( H f - -j^ ) (21) 



ferencing. We use the over-relaxed decomposition suggested by iJasatd ((1996) 
to enhance the convergence properties of the discretization of the diffusive 
term: 

^ N b - 4>c 1 , (-, d 

^— + V0. Tlf - — 

\\d\\ dMf \ d.fif 

— * 

d (d x , d y , d z ) is the vector joining the centers of the two adjacent cells (see 

Figure [2]). The average gradient V0 is interpolated from the gradients of 
these neighboring cells. The gradients of the variables at the cell centers are 
computed by a least-squares method where, for each computational cell, the 
following 3x3 linear system is solved: 

M • V0 = H (22) 

— # 

where the components of the matrix M and the vector H are written as: 

N v . N v . 

k=i \\dk\\ k= i \\dk\\ 

where N v is the number of surrounding cells. In our case, we choose to include 
all of the cells that share a vertex with cell C. Thus, in the case of a mesh 
formed by structured hexahedral cells, N v can reach 26. The contribution 
of each cell is weighted by a factor equal to the inverse of the square of the 
distance between the two cell centers „ i MO . 

The convection terms in Eq{19]are also transformed into a sum over the cell 
faces / by decomposing the surface S: 

/ pcpU.ndS = ^2( P (f>A) f Uf.n f (23) 
{ f 



3 Numerical method 



10 




where the face values cf>f require appropriate interpolation to be accurate and 
bounded. Thu s, we use the non-l inear high-resolution (HR) bounded schem e 
CUBISTA by lAlves et all u2003h in the T formulation of iNg et all (120071 ). 
where they expressed 0/ as a function of the upwind (UP) value of and its 
centered differencing (CD) value: 



lHR 



lUP 



+ T{f f D 



iUP 



(24) 



The coefficient T is determined for each face based on the local shape of the 
flow solution using the normalized variable diagram (NV D) framework and 
obser ving the convection boundedness criterion (CBC) fjGaskell and Laul . 
19881 ) . The first term of the RHS of Eq. (|24|) is accounted for implicitly, 
whil e the second term is treat ed explicitly with the deferred-correction prac- 
tice (IFerziger and Perid . Il999l). The press ure- velocity coupling is ensured by 



the SIMPLE algorithm (jPatankarl . [1980), whereas the m ass fluxes at the 



cell f aces are evaluated using the Rhie-Chow interpolation (IRhie and Chowl . 
1983T ) to avoid checker-boarding in the pressure field. Within each iteration 
of the SIMPLE algorithm, the energy equation is solved after the resolu- 
tion of the momentum equation and the Poisson equation for the pressure 
correction, and then the next SIMPLE iteration starts unless the conver- 
gence for v, p, and T is achieved. The resolution of the energy equation is 
integrated in the SIMPLE iteration to account for the high level of its cou- 
pling with the momentum equation through the body forces term and with 
the thermocapillary boundary conditions. This code is parallelized follow- 
i ng a non-ov e rlapp ing domain decomposition approach and using the Petsc 
(IBalav et all 119971 ) library. 
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4 Code validation 



This code has bee n successfully validated in some situations involving flow 



and heat transfer fjEl Omari and Le Guerl . 120101 : lEl Omari et all l201ll ) . For 
this study, we conducted additional code validations. The first supple- 
mentary validation concerns the flow in a ho rizontal channel hea ted uni- 
formly from below, which was first presented by Nicolas et al. ( 201l[ ) (steady 
Poiseuille-Rayleigh-Benard flow). Quantitative comparisons have been per- 
formed using a non-uniform mesh resolution n x = 128, n y = 176 and n z = 50. 
Figure [3] shows the longitudinal profiles of the rescaled components of the ve- 
locity and reduced temperature along the axis for the coordinates Y = 5 and 
Z = 0.5. Figure S] shows the reduced temperature field in the vertical plane 
of X = 30. As observed in Figures [3] and HI the present calculations and 
those of Nicolas et al. agree well. 

The second validation test case concerns the Marangoni convection. It is 
validated for different parameter combinations for a flow in a 2D rectangu- 
lar cavity with a free surface. The aspect ratio is fixed at T = 1, and the 
side walls are maintained at a constant temperature difference AT. The 
control parameters of the validation test ('(ISC 3X6 cLS follows: the thermo- 
capillary Reynolds number, defined as: Re t h = 7 ^ g , the Prandtl num- 
ber Pr, (the Marangoni number being Ma = Re t hPr) and the dynami- 
cal Bond number Bd = P9 ^ H , which characterizes the gravity level. This 
last parameter has the advantage of immediately giving the importance of 
the buoyancy force relative to the thermocapillary force. We perform dif- 
ferent numerical simulations for the same conditions as those described in 
Table [TJ using a uniform 121x121 grid. In Figure \5[ we plot the temper- 
ature profile obtained for Ma = 1 3 and Ma = 10 4 . The results are 

compared with those given in Ref. flZebib et all Il985h and exhibit a sat- 
isfactory c orrespondence. Other comparisons with data previously pub- 
lished by (Carpenter and Homsv . 199(J : iKuhlmann and Albensoederl . 12008 



Xu and Zebib . Il998 : Carpenter and Homsyl . 19891 ) are made: these results 
are presented in Table [TJ Our calculations are in good agreement with the 
data from the literature. 



5 Results and discussion 

To study the influence of thermal instabilities on the patterns of the flow 
and heat transfer characteristics in a horizontal channel with a free surface 
and uniformly heated from below, we examine several numerical results for 
different dimensionless parameters that control the physics of the problem. 
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Fig. 3: The longitudinal profiles of the rescaled variables Vx, Vx, Vz an d 6 (rep- 
resented by symbol s) along the axis (Y = 5, Z = 0.5) compared with the 
reference solution tNicolas et al. . 201 A) (lines). 



These parameters are the Reynolds number, the Rayleigh number, the Biot 
number and the Marangoni number. The numerical results are obtained for 
a Prandtl number fixed at Pr = 7. We are interested in using this work to 
determine the thermal instabilities that exhibit stationary longitudinal rolls 
when the Rayleigh and the Reynolds numbers are between 2500 < Ra < 
22500 and 5 < Re < 20, respectively. After a grid size-dependence study 
(see Appendix A), a mesh of resolution n x = 128, n y = 176 and n z = 50 was 
chosen for this study. 

5.1 Reynolds number effect 

As illustrated in Figure [6j the longitudinal rolls are first initiated near the side 
walls and grow gradually toward the center of the channel in the downstream 
direction. Indeed, the predominant heat transfer mechanism in these zones is 
the thermal conduction leading to a local temperature increase because the 
forced convection is weak in the vicinity of the lateral walls. This fluid heating 
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Fig. 4: The reduced te mperature field in th e vertical plane (X = 2>Q). The ref- 
erence solution {Nicolas et al . 201 A) (top) and the Tamaris code results, 
(bottom). 



induces natural convection currents near these lateral walls, causing the onset 
of the first longitudinal rolls. This onset happens near the flow inlet. The 
transverse momentum transport causes the progressive formation of other 
rolls toward the center of the channel as the fluid progress downstream. 
The rolls' dimension is nearly the same as the channel height, and the fully 
developed section contains 10 rolls. To highlight the Reynolds number's effect 
on the development of longitudinal rolls, we present the temperature field at 
the free surface for different values of the Reynolds number (5 < Re < 20 
and for Biot, Marangoni and Rayleigh numbers, respectively set to: Bi = 15, 
Ra = 5000 and Ma = 500 in Figure [7| (this value of Ma corresponds to 
weak thermocapillary effects, as will be shown in section 5.4). These flows 
correspond to a Richardson number ranging from 1 to 29. We found that an 
increased Reynolds number causes the downstream displacement of the onset 
position of the longitudinal rolls. This effect is particularly true for central 
rolls, which appear at a position X = 8 for Re = 5, while they appear at 
a position X = 32 for Re = 20. On the other hand, the rolls near the side 
walls form almost at the same position, regardless of the Reynolds number. 
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Fig. 5: Comparison of the temperature profile obtained on the free surface of a two 
dimensional square cavity for Maranqoni numb er values of Ma = 10 3 and 
Ma = 10 4 with the results of lZebib et all \l98& ) (Pr = 1). 



Figure |S] shows the variation along the channel of the transversally averaged 
Nusselt number Nu tav (cf. Eq. HBJ) for different Reynolds numbers to evaluate 
the influence of the Reynolds number on heat transfer. We note that Nut av 
first decreases and then increases to reach a maximum value before decreasing 
again slightly (behavioral characteristic of mixed convection). The maximum 
Nu tav positions correspond to the positions of the appearance of the central 
rolls. We also observe that the number of rolls (10) is not affected by Re 
in the studied range. At low Reynolds numbers, the natural convection 
improves the heat transfer for a short distance after the entry channel, causing 
Nutav to peak near X = 8. This peak has almost the same value for all Re 
numbers, but at different positions. Table [2] presents the evolution of the 
average Nusselt number Nu av over the whole heated wall with respect to the 
Reynolds number Re. We note that Nu av increases slightly with increasing 
Re, which can be explained by the transition from a mixed convection regime, 
where the heat transfer is ensured by buoyancy currents, toward a rather 
forced convection regime (relative importance of natural convection becomes 
negligible for Ri < 0.1), where the heat transfer is provided mainly by the 
forced convective flow. 
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Parameters 


Author/ Grid 


u(0,0.5) 


Nu x = .5 


Nu x =- .5 




Xu and Zebib C 1998) 


-305 


4.36 


4.36 


Re th = 10000 


Carpenter and Homsy (19901 


-296 


4.33 


4.36 


Pr = 1 


~v r 1 1 l A n l /rtnnn\ 

Kuhlmann and Albensoedcr (2008) 


-304.36 


4.364 


4.364 


Ma= 10000 


(115x115) 








Bd = 


TV 11 1 A 11 1 /riA f\ a \ 

Kunlmann and Albensoedcr (2008) 


-304.27 


4.363 


4.363 




(211x211) 
Tamaris (121x121) 


-301.99 


4.358 


4.361 


Reth = 5000 


CarDcnter and Homsv (1989) 


-179 


4.17 


4.14 


Pr = 1 


Kuhlmann and Albensoedcr (2008) 


-175.90 


4.546 


4.545 


Ma — 5000 


/ 1 1 C 1 1 r\ 

(115x115) 








Bd = 10 


Kuhlmann and Albensoedcr (2008) 


-175.97 


4.545 


4.545 




(211x211) 
Tamaris (121x121) 


-174.38 


4.543 


4.543 




Xu and Zebib (1998) 


-29.8 


6.61 


6.61 


i?e t/l = 2000 


Carpenter and Homsv (1990) 


-37.2 


6.60 


6.42 


Pr = 30 


Kuhlmann and Albensoedcr (2008) 


-24.22 


6.26 


6.26 


Afa = 60000 


(115x115) 








Bd = 


Kuhlmann and Albensoedcr (2008) 


-23.98 


6.234 


6.234 




(211x211) 
Tamaris (121x121) 


-24.64 


6.135 


6.235 



ther- 



Tab. 1: The results for the validation test case corresponding to 

mocapillary flow in a two-dimensional open square cavity com- 
pared with the results from the li teratu re Warventer and Homsi . 



199L : Kuhlmann and Albensoeden . 



Carpenter and Homsil . \1989l) . 



\2QQM; 
u(0,0.5) is the 



Xu and ZebiA , 199 1 : 



of the free surface. Nu x= q,^ and Nu x= - 
over the two vertical walls of the cavity. 



velocity at the center 
0.5 are the mean Nusselt numbers 



Re 


5 


10 


15 


20 


Ri 


28.57 


7.14 


3.17 


1.78 


Nu av 


2.31 


2.35 


2.39 


2.42 



Tab. 2: The evolution of the average Nusselt number Nu av for different Re (Bi = 
15, Ra = 5000, Ma = 500). 



5.2 Rayleigh number effect 

Figure |9] illustrates the effect on the temperature field at the free surface for 
different values of the Rayleigh number Ra. The other flow parameters are 
fixed: Bi = 15, Re = 15 and Ma = 500. One can observe that the increase 
in the Rayleigh number shifts the onset position of the longitudinal rolls 
towards the entrance of the canal. This behavior is especially true for the 
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Fig. 6: Visualization of the longitudinal rolls at selected cross sections (Re = 15, 
Ra = 5000, Bi = 15 and Ma = 500. 

central rolls, while the rolls near the side walls are formed at effectively the 
same X position (as in the previous section, the downstream displacement of 
the rolls is related to a decrease in the Richardson number). The morphology 
of the flow changes when Ra increases; in fact, the number of rolls increases, 
from 10 rolls for relatively low Ra values to 12 rolls for Ra > 12500. Figure 
mil shows the longitudinal variation of Nu tav with respect to the Rayleigh 
number. We note that Nu tav increases with increasing Ra due to the onset 
of natural convection, which appears near the entrance for large Rayleigh 
numbers (Ra > 12500), and to the increase in the kinetic energy generated 
by the existence of more longitudinal rolls. The profiles of Nu ta v also indicate 
the position where the rolls are completely developed in the cross section, 
which corresponds to the maximum Nusselt number. The slight oscillations 
observed in the Nu tav evolutions are associated with longitudinal modula- 
tions in the shape of the rolls where they start to develop. For Ra = 0, 
the case for which we do not observe roll development, the Nusselt number 
decreases continuously from the beginning of the heating zone to the channel 
outlet. 

5.3 Biot number influence 

In Figures [11] and [121 we present the transverse variation of the local Nusselt 
number at the bottom wall for various Biot numbers at X = 25 (beginning 
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Fig. 7: The reduced temperature field 9 on the free surface for different Reynolds 
numbers Re (Bi = 15, Ra = 5000, Ma = 500j. 



of the rolls' development in the cross-section) and X = 40 (all of the rolls of 
the cross-section are completely developed), respectively. The Biot number 
measures the ratio of the internal thermal resistance of the liquid film to the 
thermal resistance between the free surface and the surrounding air; a null 
Biot number corresponds to an adiabatic free surface. The other flow param- 
eters are fixed at Re = 15, Ra = 5000 and Ma = 500. We observe in these 
figures that the parietal heat transfer has a sinusoidal distribution for all of 
the studied cases. The distribution and number of peaks are related to the 
velocity distribution in the fluid generated by the ascent of hot fluid plumes 
under the action of buoyancy forces. They, therefore, reflect the organization 
of thermal instabilities into longitudinal rolls. In comparison, the results for 
a closed channel are also given (PRB flow). At X = 25 (Figure [TTT) . the 
Nusselt number evolution of the closed channel has six peaks. The two cen- 
tral peaks are of very low amplitude due to the beginning of the center rolls 
development. This is not the case when considering solutions for the free 
surface, where we see four peaks for Bi = and five peaks for Bi > 5 and 
whose amplitude is higher in the center and grows with an increasing Biot 
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Fig. 8: The longitudinal variation of the transversally averaged Nusselt number 
Nutav for different Reynolds numbers Re (Bi = 15, Ra = 5000, Ma = 
50o). 



number. When the heat exchange on the upper surface is very low (Bi ^ 0), 
the liquid film temperature increases in a generalized way, which reduces the 
thermal gradients and noticeably reduces the heat transfer at the bottom 
wall. Indeed, the convective roll number is lower for Bi = compared to 
other values. The hydrodynamic field and the Nui distribution are related; 
indeed, every peak of the Nui curve corresponds to a pair of rolls. Thus, 
8 rolls are formed for Bi = and 10 rolls form for Bi > 5. For Bi = 5, 
the two central rolls have a lower intensity than the side rolls, regardless of 
the X position, which is opposite to the behavior for Bi = 10 or Bi = 20, 
which present similar distributions with regular peaks from X = 40, but 
with slightly greater Nui values for Bi = 20 (Figure fT2l) . The presence of an 
upper wall (case of closed channel) increases the number of rolls to 12. The 
imposed temperature To at the top wall increases the thermal gradient in the 
film and should tend to improve the heat exchange. However, this favorable 
effect is counterbalanced by the viscous effects at this wall that decrease the 
fluid velocity. Thus, the Nusselt number values for the closed channel are 
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Fig. 9: The reduced temperature fields 9 on the free surface for different Rayleigh 
numbers Ra (Bi = 15, Re = 15, Ma = 500). 



lower than those obtained for Bi = 10 and Bi = 20. We also observe that 
the rolls in the case of a closed channel are regular (same wavelength), while 
in the case of the channel with free surface, the rolls close to the walls have 
a wavelength greater than those in the center. 

Figure [13] shows the longitudinal variation of the transversally average Nus- 
selt number Nu tav for different Biot numbers. For the channel with the free 
surface, the Nu tav decreases first until X = 10 and then increases gradually 
to a maximum located between X = 22 and X = 28, which corresponds 
to a position where all of the rolls are developed in the cross section. For 
the closed channel, the maximum is axially reached further downstream, at 
X = 37. Indeed, the presence of the upper rigid wall stabilizes the flow and 
delays the development of longitudinal rolls. 



5 Results and discussion 



20 



~ i r 



1 




Ra=0 
Ra=2500 
Ra= 10000 
Ra= 15000 
Ra=20000 
Ra=25000 



»oao " e s«o-©ooo-o©oo-o-©oooo-©eoo-© 



H UEH3HE 



BEH3BH&EH3HBBQB 



25 
X 



30 35 40 45 50 



Fig. 10: The longitudinal variation of the transversally averaged Nusselt number 
Nutav for different Rayleigh numbers Ra (Bi = 15, Re = 15, Ma = 500,). 



5.4 Marangoni effect 

The variation of the surface tension with the temperature (Marangoni effect 
or thermocapillary effect) generates flows at the free surface from hot areas 
towards cold areas. Figure [14] shows the effect of different Marangoni num- 
bers on the temperature fields associated with the velocity vectors at the 
cross section (X = 40) where the flow is fully developed by using different 
values of Marangoni number for Ra = 0. Setting Ra = enables us to ana- 
lyze the Marangoni effect only, without the interaction with buoyancy. We 
observe that the variation of the surface tension with temperature induces 
rolls similar to those created by natural convection. The results clearly in- 
dicate that, for Ma = 500, the flow is dominated by forced convection. The 
first longitudinal vortex rolls start to appear at Ma = 1000 where the ther- 
mocapillary effect on the basic flow becomes important. As the Marangoni 
number increases, more vortex rolls are induced and their number increases 
to reach their fully symmetrically developed state at Ma = 3000. For this 
Ma value, the maximum number of vortex rolls is 8. The rolls' dimension is 
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Fig. 11: The transverse variation of the local Nusselt number at X = 25 for dif- 
ferent Biot numbers (Re = 15, Ra = 5000, Ma = 500). 



nearly the same as the channel height. The rotation direction of the rolls due 
to the Marangoni convection is reversed relative to the rotation direction of 
the rolls due to natural convection. 

In the following section, we consider the effect of buoyancy thermal forces. 
Figure [15] shows the transverse variation at X = 30 of the transverse veloc- 
ity component V y for various values of the Marangoni number Ma and for 
Re = 15 and Ra = 5000. We observe that V y increases with increasing Ma, 
which confirms that the Marangoni effect induces a surface flow superimposed 
on the base flow. We noted also that for Ma > 1200, there is a sign inver- 
sion of the velocity that is clearly visible for the rolls near the walls, which 
is due to the Marangoni effect overcoming the dominance of the buoyancy 
thermal forces. Figure [H] presents a cross section of the flow at X = 30 and 
shows the Marangoni effect on the longitudinal rolls. Under the influences of 
thermocapillary forces, the rolls near the side walls become larger when the 
Marangoni number increases, which induces the merging of the two center 
rolls in Ma = 1200 and generates two other rolls at the sides of the channel. 
The two new rolls near the side walls then begin to grow with increasing 
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Fig. 12: The transverse variation of the local Nusselt number at X = 40 for dif- 
ferent Biot numbers (Re = 15, Ra = 5000, Ma = 500). 



Ma, causing the two center rolls to merge at Ma = 1800. While the num- 
ber of rolls decreases from 10 rolls at Ma = 1600 to 8 rolls at Ma = 1800, 
the rotational velocity of the rolls greatly increases and the rolls rotate in 
the direction imposed by the Marangoni effect. To study the influence of 
thermocapillary forces on the heat transfer, Figure [T7] shows the variation in 
the average Nusselt number at the heated bottom wall with respect to the 
Marangoni number. We note that the Nu av increases with increasing Ma 
due to the increasing velocity of the rolls under the influence of thermocapil- 
lary forces (Figure fT5|) . At Ma = 1100, we observe a decrease in the Nusselt 
number due to the beginning of the merging of the two rolls in the channel 
center. This behavior is also shown by the evolution of the transverse average 
Nusselt number along the channel (Figure US]) . Indeed, the Nusselt number 
decreases initially and then increases gradually to a maximum that indicates 
the complete occupation of the cross section by the rolls, which occurs at 
X = 24 for Ma = 800 and at X = 28 for Ma = 1100. Thus, there is a 
downstream shift of the center rolls' position with respect to the outlet of 
the channel when Ma increases between 800 and 1100. This decline contin- 
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Fig. 13: The longitudinal variation of the transverse average Nusselt number for 
different Biot numbers (Re = 15, Ra = 5000 and Ma = 500J. 



ues as Ma increases until the center rolls merge with their adjacent rolls at 
Ma = 1200 (Figure [ToT) . This reorganization of the rolls is also confirmed by 
the evolution of the velocity profile at the free surface for Ma = 1200 (Figure 
[15]) . The increase in the rolls' velocity under the action of thermocapillary 
effects promotes the heat transfer and increases the Nusselt number until a 
further decrease occurs due to the disappearance of the two center rolls at 
Ma = 1800. 



6 Conclusion 

This work presents a numerical study of three-dimensional laminar mixed 
convection within a liquid flowing in a horizontal channel heated uniformly 
from below. The upper surface is free and assumed to be flat. The side 
walls are thermally insulated. This numerical work allows us to examine the 
numerical results corresponding to a selected range of different characteristic 
non-dimensional parameters. We focused on the thermal instabilities that 
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Fig. 14: The longitudinal vortex rolls due to Marangoni effect at cross section 
X = 40. 



develop in the form of stationary longitudinal rolls. We investigated the in- 
fluence of the variation of parameters governing the physics of the problem 
on the flow patterns and heat transfer characteristics. 

The longitudinal rolls are first initiated near the side walls and grow grad- 
ually toward the center of the channel in the downstream direction. The 
rolls' dimension is nearly the same as the channel height. The Reynolds and 
the Rayleigh numbers have an important effect on the position of the rolls' 
development over the cross section. On the one hand, for a given Rayleigh 
number, the higher the Reynolds number, the greater the distance is at which 
the center rolls developed; on the other hand, for a given Reynolds number, 
higher Rayleigh numbers induce a lower distance at which the center rolls de- 
veloped. Larger Reynolds and Rayleigh numbers favorably affected the heat 
transfer by forced convection in the case of a fixed Rayleigh number and the 
development of longitudinal rolls in the case of fixed Reynolds numbers. The 
presence of an upper wall (case of closed channel) increases the number of rolls 
compared with the free surface case. The imposed temperature T at the top 
wall increases the thermal gradient in the film, which should improve the heat 
exchange. However, this favorable effect is counterbalanced by the viscous 
effects at this wall that decrease the fluid velocity and, consequently, decrease 
the heat exchange compared with the cases of Bi = 15 and Bi = 20. The 
results also show that the presence of the upper rigid wall stabilizes the flow, 
delaying the development of longitudinal rolls. Furthermore, the Marangoni 
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Fig. 15: The transverse variation at X = 30 and Z = 1 of the velocity V y for 
different Marangoni numbers Ma (Re = 15, Ra = 5000, Bi = 15J. 



effect has been investigated, and the results show that the thermocapillary 
effect on the basic flow becomes important at Ma = 1000. As the Marangoni 
number increases, more vortex rolls are induced and increase their number 
to reach a maximum of 8 rolls. The rotation direction of the rolls due to 
the Marangoni convection is reversed relative to the rotational direction of 
the rolls due to natural convection. When we combine the Marangoni effect 
with the natural convection, we observe that the thermocapillary forces start 
to overcome the dominance of the buoyancy thermal forces at Ma = 1200 
for Ra = 5000. Furthermore, increases in the Marangoni number favorably 
affect the heat transfer, but there are some drops, which corresponds to the 
disappearance of the center rolls. 



Appendix A. Mesh size-dependence study 

Several grid sizes have been tested to determine the appropriate mesh size, 
giving the best compromise between accuracy and calculation cost to ensure 
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Fig. 16: The vertical cross-section at x = 30 for different Marangoni numbers Ma 
(Ra = 5000, Re = 15 and Bi = 15J. 



that the results are grid- independent. The results are presented below for 
the case in which the control parameters are Re = 15, Ra = 5000, Bi = 15 
and Ma = 500 for the Reynolds, Rayleigh, Biot and Marangoni numbers, 
respectively. Table [3] shows the average Nusselt numbers at the bottom 
surface of the channel for increasing meshes sizes. We note that the difference 
between the results of the last two meshes does not exceed 1 7. 

Figure fT9l shows the longitudinal profile of the reduced temperature along 
the axis (Y, Z) = (5, 0.5) for the four selected meshes. We observe that 
the results remain stable and are not subjected to any influence of mesh 
at a resolution of 128x176x50. Consequently, the resolution of the mesh 
that meets the right compromise between computational time and accuracy 
chosen for this study is 128x176x50. 
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Fig. 17: The variation of the average Nusselt number at the heated wall with respect 
to the Marangoni number Ma (Ra = 5000, Re = 15 and Bi = lb). 



Meshes 


Computation time (24 processor) 


Nu av 


102x140x40 


2 h, 36min 


2.41663 


115x158x45 


4 h, 41min 


2.43877 


128x176x50 


12 h, 35min 


2.47479 


153x211x60 


51 h, 20min 


2.47639 



Tab. 3: A comparison of the average Nusselt number Nu av for various meshes 
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Fig. 18: The longitudinal variation of the transverse average Nusselt number for 
Ma = 800 and Ma = 1100. (Ra = 5000, Re = 15 and Bi = lb). 
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Fig. 19: The longitudinal profile of reduced temperature along the axis (Y, Z) = 
(5, 0.5) for different meshes. 
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